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Note: This is a working paper which will be expanded/updated frequently. All suggestions 
for improvement are welcome. The directory deleeuwpdx.net/pubfolders/pivot has a pdf 
version, the complete Rmd file with all code chunks, the bib file, and the R source code. 


1 Single Pivots 


Suppose A is an n x m matrix with a^k J 0. The k th principal single pivot or PSP transforms 
A to the n x m matrix piv fc (A) with elements 
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if i = k and j — k, 
if i = k and j J k, 
if i J k and j - k. 
if i J k and j V h- 


( 1 ) 


piv fc (A) is undefined if akk = 0 . Note that the matrix A is not necessarily square, symmetric, 
or positive semi-definite. 


If A is the matrix 


## [, 1 ] [, 2 ] [, 3 ] [, 4 ] [, 5 ] 


1 





## 

[1J 

1.00 

1.00 

1.00 

1.00 

1.00 

## 

[2,] 

1.00 

2.00 

2.00 

2.00 

2.00 

## 

[3,] 

1.00 

2.00 

3.00 

3.00 

3.00 

## 

[4,] 

1.00 

2.00 

3.00 

4.00 

4.00 

## 

[5,] 

1.00 

2.00 

3.00 

4.00 

5.00 


then pivoting on the second diagonal element, using our R function pivotOneRCO, produces 
piv 2 (A) 


## $a 


## 


[,1] 

[, 2] 

[, 3] 

[ >4] 

[, 5] 

## 

[1J 

0.5 

0.5 

0 

0 

0 

## 

[2,] 

-0.5 

0.5 

-1 

-1 

-1 

## 

[3,] 

0.0 

1.0 

1 

1 

1 

## 

[4,] 

0.0 

1.0 

1 

2 

2 

## 

[5,] 

0.0 

1.0 

1 

2 

3 


## 

## $refuse 
## [ 1 ] 0 

Note that the pivot operation is an involution, i.e. it is its own inverse, or piv fc (piv fc (.A)) = A. 

mprint(pivotOneRC(pivotOneRC (a, 2) $a , 2) $a) 


## 


[,1] 

[, 2] 

[, 3] 

[ >4] 

[, 5] 

## 

[1J 

1.00 

1.00 

1.00 

1.00 

1.00 

## 

[2,] 

1.00 

2.00 

2.00 

2.00 

2.00 

## 

[3,] 

1.00 

2.00 

3.00 

3.00 

3.00 

## 

[4,] 

1.00 

2.00 

3.00 

4.00 

4.00 

## 

[5,] 

1.00 

2.00 

3.00 

4.00 

5.00 


PSP’s are sometimes called sweeps, and the SWP() operator for symmetric matrices was 
introduced in statistical computing by Beaton (1964), section 3.2. See Goodnight (1979), 
Stewart (1998) , p. 330-333, and Lange (2010), p. 95-99, for discussions of the sweep operator 
in the context of numerical linear algebra. The version of the sweep operator that became 
popular in statistics, however, is slightly different from the one proposed by Beaton and from 
the one we use here. It is 


[swp fc (A)]y 
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a kj 

< a kk 
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if i — k and j = k , 
if i — k and j ^ k, 
if i 7 ^ k and j = k, 
if i 7 ^ k and j ^ k. 


( 2 ) 


This version preserves symmetry, which is convenient because symmetry can be used to 
minimize storage. But the symmetric sweep of $(2) is not an involution. The inverse operation 
of this sweep is the reverse sweep , which is 
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[rswp fc (-A)]jj 


'_i_ 
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if i = k and j — k, 
if i = k and j ^ k, 
if i 7 ^ k and j = k, 
if i 7 ^ k and j 7 ^ k. 


( 3 ) 


In the funcion pivotOneRCO we have a type parameter. For type equal to 1 we do a 
symmetric sweep, for type equal to 2 a reverse symmetric sweep. Thus 

mprint (pivotOneRC(a, 2, type = l)$a) 


## 


[,1] 

[,2] 

[ , 3] 

[ >4] 

[ , 5] 

## 

[1J 

0.50 

0.50 

0.00 

0.00 

0.00 

## 

[2,] 

0.50 

-0.50 

1.00 

1.00 

1.00 

## 

[3,] 

0.00 

1.00 

1.00 

1.00 

1.00 

## 

[4,] 

0.00 

1.00 

1.00 

2.00 

2.00 

## 

[5,] 

0.00 

1.00 

1.00 

2.00 

3.00 

mprint 

(pivot0neRC(pivot0neRC(a 

, 2 , type = l)$a, 2 , type = 2)$a) 

## 


[,1] 

[, 2] 

[ , 3] 

[ >4] 

[, 5] 

## 

[1J 

1.00 

1.00 

1.00 

1.00 

1.00 

## 

[2,] 

1.00 

2.00 

2.00 

2.00 

2.00 

## 

[3,] 

1.00 

2.00 

3.00 

3.00 

3.00 

## 

[4,] 

1.00 

2.00 

3.00 

4.00 

4.00 

## 

[5,] 

1.00 

2.00 

3.00 

4.00 

5.00 


Type equal to 3 is the PSP we started out with, and type equal to 4 is its transpose, the 
involution 


[q iv fc(^)]q 
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if i = k and j = k, 
if i — k and j ^ k, 
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if i 7 ^ k and j 7 ^ k. 


( 4 ) 


The default in pivotOneRCO is to set type equal to 3, and we will use this default throughout 
the rest of the paper. 


2 Sequences of pivots 

Our R function pivotRCO performs a sequence of pivots on a matrix. If both PSP’s are 
defined then piv fc (phq,(A)) = piv £ (piv A .(R)), and thus the results are independent of the 
order in which we pivot. We obtain numerical stability by always pivoting on the largest 
diagonal element in the remainder of the sequence of pivots. This also ensures the function 
gives a sensible result if the matrices involved are singular. 
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Consider the matrix we get from our previous one by setting element ( 1 , 1 ) to zero. The 
resulting matrix is still non-singular. 


## 


[,1] 

[, 2] 

[, 3 ] 

[ > 4 ] 

[, 5 ] 

## 

[1J 

0.00 

1.00 

1.00 

1.00 

1 . 

.00 

## 

[2,] 

1.00 

2.00 

2.00 

2.00 

2 . 

o 

o 

## 

[ 3 ,] 

1.00 

2.00 

3.00 

3.00 

3 . 

o 

o 

## 

[ 4 ,] 

1.00 

2.00 

3.00 

4.00 

4 . 

o 

o 

## 

[ 5 ,] 

1.00 

2.00 

3.00 

4.00 

5 , 

o 

o 

If we pivot on the first i 

four elements we get 

## 


[,1] 

[, 2] 

[, 3 ] 

[ > 4 ] 

[, 5 ] 

## 

[1J 

-2.00 

1.00 

0.00 

0.00 

0 . 

o 

o 

## 

[2,] 

1.00 

1.00 

-1.00 

0.00 

0 . 

o 

o 

## 

[ 3 ,] 

0.00 

-1.00 

2.00 

-1.00 

- 0 . 

o 

o 

## 

[ 4 ,] 

-0.00 

0.00 

-1.00 

1.00 

-1. 

.00 

## 

[ 5 ,] 

0.00 

0.00 

0.00 

1.00 

1 . 

.00 


and pivots are done in the order 4 , 2 , 1 , 3 . 

Another example is the rank 2 matrix 

mprint(a<-tcrossprod(matrix(c (1,1,1,1,1,-1,-1,1),4,2))/2) 


## 

[,1] 

[,2] 

[» 3 ] 

[ > 4 ] 

## [ 1 ,] 

1.00 

0.00 

0.00 

1.00 

## [2,] 

0.00 

1.00 

1.00 

0.00 

## [ 3 ,] 

0.00 

1.00 

1.00 

0.00 

## [ 4 ,] 

1.00 

0.00 

0.00 

1.00 

A complete pivot 

uses the order 1 , 2 , 3 , 4 and finds 

mprint (h$a) 




## 

[,1] 

[,2] 

[, 3 ] 

[ > 4 ] 

## [ 1 ,] 

1.00 

- 0.00 

0.00 

-1.00 

## [2,] 

-0.00 

1.00 

-1.00 

-0.00 

## [ 3 ,] 

0.00 

1.00 

0.00 

0.00 

## [ 4 ,] 

1.00 

0.00 

0.00 

0.00 


Note that if the function pivotOneRCO is asked to pivot on a zero diagonal element, it 
refuses and just returns its argument. Thus if zero diagonal elements are encountered, it 
is just as if the corresponding indices are not in the sequence, and we actually pivot on a 
shorter sequence. pivotRCO returns an indicator of the pivots that are skipped, which for 
our last example is 0, 0, 1, 1. 

Also note that order-invariance for pivotRCO is no longer true if pivots are skipped. Consider 

## C, 1] [,2] [,3] [ ,4] 

## [ 1 ,] 0.00 0.00 0.00 1.00 
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## 

[2,] 

0.00 

0.00 

0.00 

1.00 

## 

[3,] 

0.00 

0.00 

0.00 

1.00 

## 

[4,] 

1.00 

1.00 

1.00 

1.00 


Pivoting over all indices, but in a different order, gives different results. 

pivotRC(a, 1:4) 

## $a 


## [,1] 

[, 2] 

[, 3] 

[ >4] 

## [1,] -1 

-1 

-1 

1 

## [2,] 1 

0 

0 

0 

## [3,] 1 

0 

0 

0 

## [4,] 1 

## 

0 

0 

0 

## $pivot 
## [1] 4 1 2 

3 



## 




## $skip 
## [1] 0 0 1 

1 



pivotRC(a, 4 

: 1) 



## $a 

## [,1] 

[, 2] 

[,3] 

[ >4] 

## [1J 0 

0 

1 

0 

## [2,] 0 

0 

1 

0 

## [3,] -1 

-1 

-1 

1 

## [4,] 0 

## 

0 

1 

0 

## Spivot 
## [1] 4 3 2 

1 




## 

## $skip 

## [ 1 ] 0 0 1 1 

In the first case we compute pivr 1)4 i (^ 4 ), in the second case piv { 3 4 }(v 4 ), and these are different. 
Suppose 


with a^O. Then 


A 


a a' 

b C ’ 


piVi(A) 


~1 ^ " 

a a 

b_ _ ba 

_a a - 


( 5 ) 


5 






The Schur complement (Ouellette ( 1981 ), Zhang ( 2005 )) of ^ in piv 1 (A) is C. Since rank is 
additive over the Schur complement, we have 

rank(piv 1 (A)) = 1 + rank(C), 

and thus piv 1 (A) is non-singular whenever A is. 

For a symmetric A in ( 5 ) we have that A is positive definite if and only if a > 0 and C — — is 
positive definite (Bekker ( 1988 )). Thus if A is positive definite, all subsequent pivot elements 
are positive. If A is unit triangular, either upper or lower, then C — — is unit triangular as 
well, and thus all subsequent pivot elements are equal to one. 

Of course the choice of pivot k = 1 in ( 5 ) was for notational convenience only, the results are 
true for any k. 

3 The partial inverse 

Suppose K is a subset of N := { 1 , 2 , • • • , n} and K_ is the complimentary subset. In our 
examples we usually take, for convenience, K = { 1 , 2 , • • • , k}. The principal pivot transform 
or PPT of A on K, which generalizes the simple pivot, is 

piv x (A) = [. Ak \u , - a « kA a k J , 

\AkkA kk A kk — A kk A kk A K k 

The PPT is dehned only if A KK is non-singular. A complete PPT is dehned as piv K () for 
K = N. 

The PPT of a square matrix A is called the partial inverse by Wermuth, Wiedenbeck, and 
Cox ( 2006 ) and Wiedenbeck and Wermuth ( 2010 ). The name makes sense because if 

Akk A-k x u 

Akk Akk_ V v 

then 

Akk A kk Akk u x 

AkkA K k Akk — AkkA K kAkk y v 

And conversely. Of course all this assumes that Akk is non-singular. 

The following results are true for arbitrary index sets K and L and their complement K_ and 
L. For a proof, see Wermuth, Wiedenbeck, and Cox ( 2006 ), lemma 2 . 1 . They assume that 
all principal submatrices are invertible, which is true for positive definite matrices and unit 
triangular matrices. In that case all PPT’s are defined. 

• A = piv K (piv K (A)) 

• piv K (A) = (piv^(A))- 1 

• P iv A-(P iv A-(T>) = A -1 

• pivjpi v k (A)) = pW K (pW L (A)) 

• P^k(A) = piv A -(A _1 ) 
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• (P iv A-(^)) 1 = P iv x(^ x ) 

Wermuth, Wiedenbeck, and Cox ( 2006 ) summarize the main properties of partial inversion in 
their theorem 2 . 2 . Suppose K, L, and M partition the n indices of a square matrix, and the 
PPT’s are defined Then 


• Commutativity: piv A -(piv L (A)) = piv A -(piv L (A)) = piv // (A). 

• Exchangeability: [pi v K (A)] KuLtKuL = pi v K (A KuL>KuL ). 

• Symmetric Difference: piv Aui (A)piv LuM (A) = piv KuM (A). 

If A is not square, then some of these results still make sense. For both a broad n x m matrix 
A, with n < m, and a tall A, with n > m, we can pivot on the diagonal elements and the 
defining equations for the principal pivot transform still apply. 


Note that a PPT can be computed by a sequence of PSP’s, provided all PSP’s on the way 
are defined. But consider the example 


## 


[,1] 

[, 2] 

[, 3 ] 

[, 4 ] 

## 

[1J 

0.00 

1.00 

1.00 

0.00 

## 

[2,] 

1.00 

0.00 

0.00 

1.00 

## 

[ 3 ,] 

1.00 

0.00 

1.00 

0.00 

## 

[ 4 ,] 

0.00 

1.00 

0.00 

1.00 


Now piv 1 (A) and piv 2 (A) are both not defined, and thus our function pivotRC(a, 
does nothing, while 



' 0 

1 

0 

-r 

P iv {l,2}(4> = 

1 

0 

-1 

0 

0 

1 

1 

-1 


1 

0 

-1 

1 


1 : 2 ) 


4 Inverse and determinant 


For a square matrix 


A 


a a' 
b C ’ 


with a ^ 0 , by Schur’s determinant formula, 


det(A) = a det(C-), 

a 


and thus det(A) is the product of the pivots in a complete PPT. 

If A is square and non-singular, then carrying out a complete PPT produces the inverse of 
the matrix, or piv Ar (A) = A -1 . This assumes, of course, that a complete PPT can actually 
be carried out and we do not run into a zero pivot. 

If A is singular, say of rank r, then we can permute rows and columns so that we have the 
leading principal submatrix of order r is non-singular. This easily follows from the pivoted 
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LU decomposition, see Stewart ( 1998 ), theorem 2 . 13 . If 


A 


B C 


D E 


with B an r x r matrix such that rank(R) = rank(A) = r, and a pivot of A over the first r 
elements is possible, then 




' B~ l —B~ l C 
DB~ l 0 


This implies that A (piv{ lj2) ... )T .j(-A)) A = A, i.e. piv{ 12; ... jr }(A) satisfies the first Penrose 
condition, i.e. is a (l)-inverse of A. 


5 Least Squares 

Statisticians like the piv() and swp() operators, because of their connection with linear least 
squares regression, and with stepwise regression in particular (Jennrich ( 1977 )). In addition 
there are connections with covariance selection (Dempster ( 1972 )) and with graphical models 
(Wermuth, Wiedenbeck, and Cox ( 2006 )). 

For the least squares problem of minimizing (y — Xb)'(y — Xb ) the normal equations are 
X'Xb = X'y. Suppose X = (Ad | X 2 ), with rn i and m2 columns, and consider the matrix 

X(Ad X[X 2 X\y 
A = X'Ad X'X 2 X' 2 y 
_ y'X\ y'X 2 y'y 

Pivoting on the m\ elements in X[X\, assuming it is invertible, gives 

(X(Ad)” 1 -{X[X\)~ l X[y 

Piv {1)2 ,., m} (^) = X' 2 Xi{X[Xi)~ l X'(/-Ad( X'yX^X'jXi X' (/ - X^X'X^X'Jy . 

y'XAX'Xt )- 1 y\I - X^X'yX^-'XtiX* y\I - X 1 (X l 1 X 1 )~ 1 X l 1 )y 

All submatrices in the PPT give statistics which are relevant in the context of the least- 
squares problem. We have the regression coefficients, their dispersion matrix, the residual 
sum of squares, and the partial covariances of A" 2 given X\. Moreover an additional PSP on 
one of the m 2 elements from X 2 will add a variable to the regression, and repeating a PSP 
from the first mi will remove that variable from the regression. ^Code 

The basic PSP routine pivotOneCQ is written in C. If a pivot is zero, it is skipped. Performing 
a number of pivots to get a PPT is done in pivotCO, also in C. The ,C() interface is used 
in R to write wrappers pivotOneRCO and pivotRCO for the C routines. Note that the ggm 
package from CRAN (Marchetti, Drton, and Sadeghi ( 2015 )) has a swpO function that uses 
solve () from base R to compute the partial inverse. Pivoting on principal submatrices to 
compute the partial inverse will generally be faster than building up the partial inverse from 
pivoting on a sequence of diagonal elements like we do. 
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In previous work ((???)) we implemented sweeping by writing a simple R wrapper for the 
fortran routine in AS 178 (Clarke ( 1982 )). We improve this by switching to principal pivoting, 
to C, and to incorporating reordering the remaining pivots for numerical stability, similar to 
the suggestions of Ridout and Cobby ( 1989 ). 

5.1 R Code 

dyn.load ("pivot.so") 


pivotOneRC <- function (a, k, type = 3, eps = le-10) { 
n <- nrow (a) 
m <- ncol (a) 
h <- 

• C( 

"pivotOneC" , 
a = as.double (a), 
k = as.integer (k), 
n = as.integer (n), 
m = as.integer (m) , 
type = as.integer (type), 
refuse = as.integer (0), 
eos = as.double (eps) 

) 

return (list(a = matrix(h$a, n, m), refuse = h$refuse)) 

> 

pivotRC <- function (a, ind, type = 3, 
n <- nrow (a) 
m <- ncol (a) 
p <- length (ind) 
done <- rep (0, p) 
h <- 

• C( 

"pivotC" , 

a = as.double (a), 
ind = as.integer (ind), 
jnd = as.integer (rep (0, p)), 
done = as.integer (rep (0, p)), 
n = as.integer (n), 
m = as.integer (m), 
p = as.integer (p), 
type = as.integer (type), 
skip = as.integer (rep (0, p)), 
eps = as.double (eps) 


eps = le-10) { 
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) 

return (list(a = matrix(h$a, n, m), pivot = h$jnd, skip = h$skip)) 

> 

5.2 C Code 

void pivotOneC (double *, int *, int *, int *, int *, int *, double *eps); 

void pivotC (double *, int *, int *, int *, int *, int *, int *, int *, int 

void mprint (double *, int *, int *); 

void pivotOneC (double *a, int *kk, int *nn, int *mm, int *type, int ^refuse 

int i, j, ik, k j , i j , kp, sOO, sOl, slO, n = *nn, m = *mm, k = *kk, tp = 

double pv, ee = *eps; 

^refuse = 0; 
if (tp == 1) { 
sOO = -1; 
sOl = 1; 
slO = 1; 

} 

if (tp == 2) { 
sOO = -1; 
sOl = -1; 
slO = -1; 

> 

if (tp == 3 ) { 
sOO = 1; 
sOl = -1; 
slO = 1; 

> 

if (tp == 4 ) { 
sOO = 1; 
sOl = 1; 
slO = -1; 

> 

kp = MINDEX(k, k, n); 

pv = a [kp] ; 

if (fabs (pv) < ee) { 

^refuse = 1; 
return; 

> 

for (i = 1; i <= n; i++) { 
if (i == k) continue; 
ik = MINDEX (i, k , n); 


, double *e 

double *ep 
*type; 
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for (j = 1; j <= m; j++) { 
if (j == k) continue; 
kj = MINDEX (k, j, n); 
ij = MINDEX (i, j, n); 
a[ij] = a[ij] - a[ik] * a[kj] / pv; 

} 

> 

for (i = 1; i <= n; i++) { 
if (i == k) continue; 
ik = MINDEX (i, k, n); 
a[ik] = slO * a[ik] / pv; 

> 

for (j = 1; j <= m; j++) { 
if (j == k) continue; 
kj = MINDEX (k, j, n); 
a[kj] = sOl * a[kj] / pv; 

> 

a[kp] = sOO / pv; 


void pivotC (double *a, int *ind, int *jnd, int *done, int *nn, int *mm, int *pp, int *t 
int ipiv, kpiv, lpiv, refuse; 
int n = *nn, m = *mm, P = *pp; 
double fmax, fpiv; 
for (int 1=0; 1 < p; 1++) { 
jnd[l] = 0; 
done[1] = 0; 
skip[l] = 0; 

> 

for (int 1 = 1; 1 <= p; 1++) { 
fmax = -1.0; 

for (int k = 1; k <= p; k++) { 
kpiv = ind[k - 1] ; 
if (done[k - 1] == 1) continue; 
fpiv = fabs (a[MINDEX (kpiv, kpiv, n)]); 
if (fpiv > fmax) { 
ipiv = k; 
lpiv = kpiv; 
fmax = fpiv; 

} 

} 

done[ipiv - 1] =1; 
jnd [1 - 1] = lpiv; 
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pivotOneC(a, felpiv, &n, &m, type, ferefuse, eps); 
skip[l - 1] = refuse; 

> 

> 


void mprint (double * a, int *nn, int *mm) { 
int i, j, ij, n = *nn, m = *mm; 
for (i = 1; i <= n; i++) { 

for (j = 1; j <= m; j++) { 
ij = MINDEX (i, j, n); 
printf (" %5.3f ", a[ij]); 

} 

printf ("\n"); 

> 

> 

6 NEWS 
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